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Abstract 

We calculate thermal production of u, d, s, c and b quarks in ultra-relativistic heavy 
ion collisions. The following processes are taken into account: thermal gluon decay 
{g — > u), gluon fusion {gg —>■ if), and quark- ant iquark annihilation (jj ii), where i and 
j represent quark species. We use the thermal quark masses, rnfiT) m] + (2(7^/9)T^, 
in all the rates. At small mass {mi{T) < 2T), the production is largely dominated by 
the thermal gluon decay channel. We obtain numerical and analytic solutions of one- 
dimensional hydrodynamic expansion of an initially pure glue plasma. Our results show 
that even in a quite optimistic scenario, all quarks are far from chemical equilibrium 
throughout the expansion. Thermal production of light quarks {u, d and s) is nearly 
independent of species. Heavy quark (c and b) production is quite independent of the 
transition temperature and could serve as a very good probe of the initial temperature. 
Thermal quark production measurements could also be used to determine the gluon 
damping rate, or equivalently the magnetic mass. 
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1. Introduction 



Recently, there has been much interest in calculating the production rates for both 
massless and massive quarks in ultra-relativistic nuclear collisions Unfortunately, 
there have been no rigorous thermal field theory calculations to compare to approximate 
results, probably because the "naive" production rates diverge for massless quarks. The 
production rates can be regularized by summing over hard thermal loops 0, § , following 
the treatment of Braaten and Pisarski 0. This is a difficult calculation, and has not 
been done so far. 

The closest approach to a full calculation has been to use the thermal quark masses 
for massless quarks to regulate the divergent matrix elements However, the thermal 
corrections for massive quarks are not included, and these can be important when m ~ 
T. We also take into account the thermal gluon decay process |10[, which was never 
considered in previous works. In this way, we can see if the quark production rates in 
non-thermal QCD simulations of collisions are similar to those in thermally-equilibrated 
quark-gluon matter at equivalent temperatures. 

In this paper, we calculate thermal quark production from thermally-equilibrated 
quark-gluon matter in ultra-relativistic nuclear collisions. Our initial state is a high- 
temperature gluon plasma (GP), with gluons in thermal and chemical equilibrium but 
no quarks or antiquarks. Produced quarks are assumed to be in thermal equilibrium, as 
thermal equilibration times are short, but not in chemical equilibrium. Production of 
low-mass (mj < 3.5T) quarks is dominated by thermal gluon decay {g — >• iz), while gluon 
fusion {gg ii) is the most important production mechanism for high mass (m^ > 3.5T) 
quarks. We also include quark-antiquark annihilation (j j it) , where i and j represent 
different quark species, although this process is shown to be unimportant to our analysis. 

In section we calculate the rates for the processes that we include in our simulation. 
In section ^, we describe our simulation and calculate effects of the rates on quark 
chemical equilibration and thermal quark production. We summarize our results in 
section 0. 



2. Quark production rates 

In this section we calculate the rates for the various processes. We begin with the 
production rate from quark-antiquark collisions, Rjj^n. We then calculate the rates 
from gluon fusion and thermal gluon decay, Rgg^n and Rg^u respectively. 

For most of the collision, the quark densities are below their equilibrium values, so the 
quark chemical potentials fii, fij < 0. In the absence of baryons, the antiquark densities 
are equal to the quark densities, so we also have fj,j = fii < 0. As the initial densities are 
near zero, the initial values —fii, —Hj ^ T, in which limit the rate Rjj^ii can nearly be 
expressed in closed form. The exact rate near equilibrium is not so important for our 
results, so we use these far-from-equihbrium rates throughout this paper. 

There does not yet exist any consistent formalism which is able to deal with systems 
out of chemical equilibrium. In particular, merely substituting non-equilibrium quark 
and antiquark distributions for equilibrium Fermi-Dirac distributions is certainly flawed 
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beyond the one-loop level. Indeed, in chemical equilibrium, the quark chemical potential 
is opposite to the antiquark chemical potential. This guarantees the micro-reversibility 
relations, which are at the basis of a reliable perturbative expansion |Tl| . In the present 
work, we use non-equilibrium densities (with the same chemical potential for quarks and 
antiquarks) at the tree-level only. 

Even though quarks are not in chemical equilibrium, they thermalize very fast with 
the gluons. The gluon correction to the thermal quark mass is 



2 I 2(7'^ 2 
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in the high-frequency limit, which we need for most of our calculations. As there are 
very few quarks in the plasma, we can simply drop the thermal quark contribution to 
the thermal quark mass. This amounts to using the same approximation as in the rate 
calculations, — yUj, —fij ^ T. 

As ~ 2, this gives us minimum masses of approximately T. These minimum masses 
are large enough that we can use the high-mass limits of quark production rates in our 
simulation. 



2.1. Quark-antiquark collisions 

We begin by calculating the rate for the process jj 
d^Pj d^pj d^pi d^pi 
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where \M.^ is summed over initial and final spins and colors, and 

/fd(^0 = /fd(^^) ^ 

In the limit — yUj, — yUj T, eq. (H) reduces to 

,2^1 j/T f d.^Pj d^Pj 
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Following Ref. M, we evaluate all but 5 of the integrals trivially, obtaining 
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where is evaluated with 
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Here M{m) and -Po(Po) are the mass and half the energy difference of the heavier (hghter) 
quark and antiquark, while u and q are the pair energy and momentum in the matter 
rest frame. The limits of integration, 
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arise from the kinematic constraints. 

For the process jj it, the summed squared matrix element is 

(M^ + - t)2 + (M^ + - + 2{M^ + m^)^ 
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Inserting the summed squared matrix element (|TT]) in the rate (^ and evaluating the 
integrals over 0, po and Pq; we obtain 
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In general, the remaining two integrals must be evaluated numerically. The rate can be 
evaluated analytically when M <^ T, where 
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The two limits are almost the same if evaluated with a mass M ~ T; the ratio of the 
two rates is 3itM/AT if the exponential is omitted and the lighter mass is set to zero. We 
use only the result for the heavy quark limit in our simulation, as the difference in the 
rates for the light quarks is then not very large if we give them thermal masses. Also, 
the total rate in this low-mass region is dominated by the gluon decay process. 

In the limit of low quark density, the forward and reverse rates differ only in the 
chemical potential for the initial-state quarks, so we find 
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We use the convention throughout this paper that 

SRx^Y = Rx^Y — Ry^x- 
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2.2. Gluon fusion 



The rate Rgg^ij has not yet been properly calculated in perturbative thermal QCD. 
Naive calculations diverge for massless quarks, so the rates must be regulated by includ- 
ing hard thermal loops. This has not yet been done; only the dominant contribution for 
small g has been calculated [ IC ] . We again approximate the full calculation by using the 
low-density thermal quark masses to regulate the production rates. 

The rate for the process gg — n is Q 

_ 1 f d^pi d^p2 d^Pi d^pi -4 4, _ ^_ 

^99^n - 2 j (27r)32Ei (271)32^2 {2Tif2E, {2ny2E^ ' ^ 

xEI-^..-//be(^i)/be(^2), (17) 
where the subscripts 1 and 2 refer to the initial-state gluons and 

/be(^) = ^EJ^Y (18) 

Here we have assumed that the heavy quark density is zero, so that there is no Pauli- 
blocking of final states. We proceed in the same manner as in our calculation of Rj^^ij, 
obtaining 

1 POO /.(a;2_4^2^1/2 

Rgg^iJ = .^.o ^6 duj dq dpo / dp^ 
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where u, q and are defined as before, 

po = E1-E2, (20) 

p', = E,-E,, (21) 

and the limits of integration are 

P* = |, (22) 

A - 1(1-^) ■ P3) 

Up to this point, we have made no approximations (except for ignoring the final- 
state Pauli blocking). Now, we want to perform the calculation analytically for the case 
Mi ^ T. In this case, we again define uj^ = u — 2Mj, obtaining 

I roo /■(4A/ia;.+(^2)l/2 
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It is clear from eq. (||) that cu^ ~ T, so that Po ^ q ^ (AMiUj^y/'^ ~ {MiTf/'^ < M^. 
Thus, the argument of f-Q^j is always large, so we can replace the Bose distributions by 
exponentials. Similarly, we can show that 
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t = u 



4M^{1 + 0[T/Mi]}, 
-M^{l + 0[{T/M^y/^]}. 



We can then easily evaluate the integrals to obtain 
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where M. is evaluated at s = AMf,t = u = — M?. Note that the fractional correction in 
R is of higher order in T /Mi than the fractional correction in t or m; this is because M. 
is invariant under the transformation t ^ u, and because s + t + u = 2Mf the correction 
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Inserting this into eq. (P7[), we find 
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Again, the forward and reverse reactions differ only in chemical potential, so we obtain 
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2.3. Thermal gluon decay 

In a previous work |10|, we found that the thermal gluon decay into quark-antiquark 
pairs is the leading contribution in the perturbative expansion, at least for massless 
quarks. As it is also relevant for massive quarks, we include it in our numerical simula- 
tion. 

A diagrammatic analysis reveals that one of the diagrams considered previously in 
the gluon fusion process, namely the s-channel, is in fact included in the thermal gluon 
decay set of diagrams. There is no over-counting, though, as we do not integrate the 
same phase space region in the two cases. 

We evaluate the integrals for Rg^ii in the limit Mi/T ^ l,g. We begin with the 
expression for the quark production rate from transverse gluon decay. 
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where 7^ is the gluon damping rate and IIt{Q) is the transverse projection of the po- 
larization tensor (the square of the thermal mass for transverse gluons). In the region 
where the integrand is large, 

7g = lasT In — + 0{asT) and Renr(Q) = ^oo^o = '^t^c^sT^- (32) 
2 as 2. 

In the high-mass and small coupling limit, eq. (|3lD reduces to 
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where ml = tu^ — g^. The integral falls off exponentially fast with increasing u, so we 
expand in u^: = u — 2Mi and 6ml = ^l ~ 4:M^: 
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times [1 + OiT /2Mi)\. The IIt and 7^ in the denominator provide correction factors 



[1 -|- nT/2Mf] and [1 — 7g/Mf ] respectively, and these can also be dropped if Mj is large 



enough, giving 
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where z = 5ml/ AMiUj^,. These integrals can be evaluated analytically, giving 
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The longitudinal case is the same except for a factor of 1/2: 
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The total net production rate is 
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This rate is significantly enhanced for light quarks, however, so we use the exact rate in 
our simulation. 
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2.4. Comparison of the rates 



We can now numerically compare the different rates that we have obtained and com- 
ment on our different approximations. In Fig. 1, we show the quark: production rate as a 
function of the quark mass m, taking /i — — oo for outgoing quarks. The thermal gluon 
decay dominates for m < 2 T and is about ten times larger than Rjj^n and Rgg^ij, at 
m ~ 0.5 T. For the latter two processes, we have used the high mass formulas, eqs. (0) 
(with yUj = 0), and ( ]29|) respectively. As rrii 0, the quark-antiquark annihilation rate 
tends asymptotically to 4/(37r^)a| ~ 3 x 10~^, according to eq. (p!3[) . Therefore, it is 
clear that we can reliably use the high- mass formulas down to m ~ 0.5 T. The ther- 
mal masses being of order T, we see that we do not need in our simulation anything 
other than those high-mass formulas. It makes irrelevant the only uncertainty in our 
calculation regarding the process Rgg^n, which has not been calculated rigorously in the 
low-mass region. Note also that we have neglected the thermal gluon mass in calculating 
Rgg^ii, which is justified as this process starts to dominate only when m > 3.5 T, which 
is typically large compared to the thermal mass. 

For the thermal gluon decay, the asymptotic expressions (|36D and (|37|) turn out to be 
close to the exact value only when m > 2 T. We have not been able to make a good fit 
of the numerical curve for any values of g. We therefore use the exact expression for the 
gluon decay, eq. (^9]), in our simulation. 



3. Quark production in a nuclear collision 

Our simulation works as follows. We assume that at initial proper time tq the glu- 
ons are in thermal and chemical equilibrium at initial temperature To but no quarks 
are present. We then perform boost-invariant one-dimensional (no transverse compo- 
nent) hydrodynamic expansion, conserving all quark numbers but not gluons, then allow 
number-changing processes explicitly via chemical reactions. All particles are assumed 
to be in thermal equilibrium (and the gluons in chemical equilibrium) throughout the 
collision. 

The evolution equations are 

de e + P 



dr T 
drii —Hi 

dr T 



(40) 



+ SRg^U + SRgg^U + ^Rjj^it- (41) 



The first terms of eqs. (^0]) and (|4T|) give the hydrodynamic evolution, and the remaining 
three terms of eq. (^Tj) give the chemical evolution of the quarks. 

The temperature and chemical potentials are recalculated after each time step from 
the densities and chemical potentials, using the thermal equilibrium conditions 

3 r°° 

n,{T,fii) = n,{T,fii) = — dk f^jy {Ei{k)) , (42) 
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The first contributions to the energy density and the pressure are from the gluons. All 
of the above integrals are evaluated numerically to better than 1% accuracy. 
We use the same running coupling constant as in Ref. f^: 

"5 = — r^,.„T (45) 



(33 - 2n/) In [QVA 

Here ^/ is the number of light flavors for the given reaction (4 for h production, 3 for 
all other reactions), Q"^ is the mean momentum transfer, and A is the renormalization 
scale. We use A = 300 MeV for our calculations, instead of 400 MeV as in Ref. 0; this 
gives larger values of and thus higher reaction rates. We take the ansatz 

Q2 = m\T) + 9T2, (46) 

which is correct in the high-mass limit and for a Boltzmann distribution in the low- 
mass limit. We have checked that this as-behavior is in good agreement with recent 



calculations on the running coupling constant at finite temperature |jT3|. It gives a 
rather weak A-dependence. 

For our first simulation, we use conditions for an Au+Au reaction at = 200 
GeV/nucleon: 

To = 0.063 fm/c, (47) 

To = 950 MeV, (48) 

The initial time is somewhat arbitrary; we picked these numbers because they give 
the same mass spectrum for thermal dilepton production as a more detailed collision 



simulation |T^. For quark masses we use = nid = 0, m<( = 0.2, rric = 1.5 and mi, = 5 
GeV. We run until r = 11 fm/c, corresponding to T ^ 150 MeV. We refer to this as our 
optimistic scenario because the extremely short equilibration time (and correspondingly 
high temperature) result in large values of heavy quark production. 

In Fig. 2 we show the ratio of density to chemical equilibrium (/i = 0) density for each 
quark species. This ratio gives a good measure of the degree of chemical equilibration 
for the light quarks (m, d and s), as they are produced throughout the evolution. Note 
that, because of the thermal mass effects, thermal production of light quarks is nearly 
independent of species. Heavy quarks (c and h) are produced only in the early stages, 
when their densities are far below equilibrium values. In the later stages their densities 
actually exceed the equilibrium values, but at this point they are no longer interacting 
chemically so they are not equilibrated (they pass the equilibrium values because the 
equilibrium densities fall exponentially fast with increasing r). 

We see from Fig. 2 that all quarks are far from chemical equilibrium throughout the 
collision, in agreement with the conclusions of Geiger and Kapusta |]15[. However, our 



results differ somewhat from theirs. The equilibration times are longer, so that ratios 
of density to equilibrium density rise more slowly at first. Also, the ratios of density to 
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equilibrium density are always increasing throughout our simulation; this is because all 
quark densities are less than their equilibrium values except for the heavy quarks after 
freezeout. 

In Figs. 3 and 4, we show the quark (or antiquark) production per unit rapidity at 
y = for each species as a function of Tc, using respectively our optimistic scenario and 
a more standard scenario with the same entropy density: 

To = 0.2fm/c, (49) 
To = 650 MeV. (50) 

The rapidity density is simply related to the density at Tc'. 

dNi/dy\y=Q = ^™ilT=re ■ (51) 

where A is the cross-sectional area of the plasma. We use A — 150 fm^ for a central 
Au+Au coUision. 

These results can be understood analytically. For M ^ Tq, the heavy quark approxi- 
mations arc clearly justified. The dominant production mechanism is then gluon fusion, 
as Rgg^Qql Rg^QQ oc M/T. Quark collisions are unimportant, as all quark densities are 
small throughout the collision. 

Because the quark densities are small, the thermal behaviour is dominated by gluons, 
so that to a good approximation 

tT^ = ToTl (52) 
as expected from entropy conservation. We can then estimate the quark densities: 



-(r) = 1 [ dr' r' Rgg^,,{r') = ^^Mr^T.^T^ j^' ^e^-Z^'. (53) 



We then obtain the contribution to the rapidity density from thermal gluon-gluon colli- 
sions in the plasma. 
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Wc can make a similar calculation of the contribution from thermal gluon decay, 
obtaining 
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When 2M/Tc — 2M/Tq » 1, the heavy quark production in the plasma is independent 
of T.: 
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(56) 



This is because heavy quark production occurs only at the earhest (and hottest) part of 
the collision. 
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In Fig. 5, we compare c and b quark production in our simulation with the analytic 
formula for heavy quarks (pUI). We calculate as from eq. ( ^Sj) with T = Tq, and use 
this value to calculate 7^ and the thermal mass; this works well because most of the 
heavy quark production occurs near Tq. The analytic result is a bit high for the b quarks 
because the quark production lowers the temperature faster than for a pure glue system. 

The situation for massless quarks is somewhat different. Here the dominant production 
mechanism is decay of thermal gluons. We thus find 

Rg^,^ = KT\ (57) 



.(58) 



K = / / dujujfBE{oj) / ds [s + 2i 

'( 2_\( _^\X'^ \ I 1 

Here v is the ratio of the thermal quark mass to T: 



\ 1/2 

(59) 



We then obtain 



Similarly, we find 



2' 



7 f T 



2' 



^^f /^^ = ^2^'s''^~"'^^oToX-' - • (61) 

Note that, for massless quarks, when Tq ^ Tc we find that 

dNjdy = (^^K + ^c^l^e-'") Ar^T^T^-'. (62) 

Light quark production in the plasma is not independent of Tc. The difference between 
analytic and simulated results can be seen in Fig. 6. The agreement is not so good as in 
the heavy quark case but still remains within a factor of two. The other main difference 
from the heavy quark production is that the best value of as to use is not given by 
the value at Tq, as the light quark production is not dominated so much by the early 
stages. Of course, this quantity is also relatively meaningless for the u and d quarks, 
as these are copiously produced in normal hadronic (confined) matter at temperatures 
near T^. However, the s quarks are not so easily produced in confined matter, so their 
production in the plasma is of interest. Also, the analytic result serves as a useful check 
of the simulation. 

Let us finally list some remarks: 

- The sensitivity to the choice of the renormalization scale A is rather weak. Our results 
vary less than 30% for A = 300 — 400 MeV At the beginning of the expansion, 05 is 
always small and perturbative QCD applies, so our calculation of heavy quark production 
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is reliable. Light quarks are still being produced at the end of the plasma phase (near 
Tc), where large values of as can lead to significant quark-antiquark production and 
annihilation. The subsequent hadronization also introduces large uncertainties in light 
quark production, especially for the u and d fiavors. 

- An exciting feature is that the gluon decay process is responsible for much of the 
heavy quark production. The simulations show that, with Tq = 950 MeV, it accounts 
for 70% of charm production and 50% of bottom production. If Tq = 650 MeV, it still 
represents 60% of charms and 40% of bottoms which are produced during the expansion. 
These numbers agree well with the analytic estimate (0). Therefore, it seems very 
realistic to measure the gluon damping rate, or equivalently the magnetic mass, at future 
heavy-ion colliders if a plasma with sufficiently high temperature can be formed. 

4. Conclusions 

We have calculated thermal quark production in ultra-relativistic heavy ion collisions. 
The following processes have been taken into account: thermal gluon decay {g it), 
gluon fusion {gg — » it), and quark-antiquark annihilation {jj ii). We have used the 
thermal quark masses in all the rates. Quark production is dominated by the gluon 
decay channel for small quark mass, and by gluon fusion for high mass. 

We have performed a one-dimensional hydrodynamic expansion of a plasma, initially 
free of quarks. Our results show that even in a quite optimistic scenario, that is for a 
simulation with initial time tq = 0.063 fm/c and an initial temperature To = 950 MeV, 
the light quarks {u, d and s) do not reach equilibrium at the end of the expansion. 
Furthermore, due to thermal mass effects, light quark (m, d and s) production is nearly 
independent of species. Heavy quark (c and h) production is quite independent of the 
transition temperature and could serve as a very good probe of the initial temperature. 
Thermal quark production could also be used to determine the gluon damping rate, 
provided the initial heavy quark density created in parton collisions (before the formation 
of the plasma) is not too large. This is equivalent to determining the magnetic mass. 
All these numerical results have been checked to be in very good agreement with simple 
analytic estimates, at least for the heavy quarks. For light quarks, the agreement between 
analytic and numerical results is within a factor of two. 

Obviously, these results should have deep implications for other signals, as thermal 
photon and dilepton emission, for which quark chemical equilibrium is usually assumed 
|TT|,|TB|. At the moment, it is not clear whether the relatively high initial temperature 
can counterbalance the effect of having four to five times less quarks than in a complete 
equilibrium picture. 

Finally, the large rates of heavy quark production seen in the parton cascade model 
do not appear in our simulation. This indicates that they are probably a result of 
high-energy parton collisions, and not from any quark-gluon plasma that is included 
(implicitly) in the parton cascade model. 
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Figure captions 



1. Comparison of quark production rates in tlie different channels. We use the com- 
plete result for the gluon decay, eq. (|39D, and high-mass formulas for gluon fusion 
and quark-antiquark annihilation, eqs. ( p9D and ( |T^ ) respectively, with as = 0.32. 
Outgoing quarks are unthermalized, and incoming quarks {u and d) are in thermal 
equihbrium = /x^ = 0). 

2. Degree of chemical equilibration in a central Au+Au collision with ^/s = 200 
GeV/nucleon for the optimistic scenario (tq = 0.063 fm/c, Tq = 950 MeV). 

3. Rapidity densities in a central Au+Au collision with ^/s = 200 GeV/nucleon for the 
optimistic scenario (tq = 0.063 fm/c, Tq = 950 MeV). 

4. Rapidity densities in a central Au+Au collision with ^/s = 200 GeV/nucleon for the 
standard scenario (tq = 0.2 fm/c, Tq = 650 MeV). 

5. Comparison of simulation with analytic results for heavy quark production, eq. ([56|), 
for the optimistic scenario (tq = 0.063 fm/c, Tq = 950 MeV). 

6. Comparison of simulation with analytic results for light quark production, eq. ([HT|), 
for the optimistic scenario (tq = 0.063 fm/c, Tq = 950 MeV). 
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